A New Algorithm for Image Reconstruction and Image Noise 
Analysis in Computed Tomography 

Field of the Invention 

[00011 The field of the invention relates to computed tomograpy and more 
particularly to methods of reconstructing images. 

Background of the Invention 

[0002] The U.S. Government has a paid-up license in this invention and the right in 
limited circumstances to require the patent owner to license others on reasonable terms as 
provided for by the terms of Grant No. EB00225 and CA70449 awarded by the National 
Institutes of Health. 

[0003J In computed tomography (CT), the fan-beam filtered backprojection (FFBP) 
algorithm is typically used for reconstruction of images directly fi-om fan-beam data. In single- 
and multi-slice helical CT, the FFBP algorithm may also be used for reconstructing multiple 
slices of two-dimensional (2D) images firom fan-beam data that are converted from measured 
helical data. The FDK algorithm may be used for cone-beam data. It has been known that the 
spatially variant weighting factor in the FFBP algorithm can significantly amplify data noise and 
sample aliasing in situations where the focal lengths are comparable to or smaller than the size of 
the field of view (FOV). One previously developed hybrid algorithm reconstructs images &om 
parallel-beam data that are converted from the acqiured fan-beam data. Because the previous 
hybrid algorithm involves no weighting factor similar to that in the FFBP algorithm, it is 
generally less susceptible to data noise and sample aliasing and is computationally more efficient 
than is the FFBP algorithm. However, because the previous hybrid algorithm invokes an explicit 
one-dimensional (ID) interpolation for converting fan-beam projections to parallel-beam 
projections along the radial direction, it can lead to reduced unage resolution when fan-beam 
samples along the radial direction are sparse. 

[0004] A new hybrid algorithm is described herein that retauis the favorable 
resolution property of the FFBP algorithm and the favorable noise property of the previous 
hybrid algorithm while eliminating their shortcomings. Analytic formulas are also derived for 
evaluating variances in unages reconstructed by use of the FFBP, the previous hybrid, and the 
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new hybrid algorithms in their discrete forms. Such theoretical formulas not only may provide 
insights into the algorithms' precision in estimation tasks, but also may be used for assessing 
their performance by use of mathematical model observers in detection/classification tasks. 
Computer-simulation studies for quantitative evaluation of resolution and noise properties of 
these algorithms are provided. Using moduli of Fourier transforms of reconstructed images, 
image-resolution properties may be compared and possible distortions identified. From a large 
number of reconstructed noisy images, empirical image variances are computed and compared 
with the derived theoretical formulas. Numerical results of these simulation studies confirm that 
the proposed new hybrid algorithm combines the favorable resolution property of the FFBP 
algorithm and the noise property of the previous hybrid algorithm while eliminating their 
shortcomings. Empirical image variances also validate the novel theoretical formulas for image 
variances. The analysis contained herein has also be extended to half-scan CT, cone beam CT, 
asynmietric scanning and high resolution scanning techniques. 

Smnmary 

[0005] A method and apparatus are provided for reconstructing a tomographic image 
from fan-beam or cone-beam data. The method includes the steps of collecting fan-beam or 
cone-beam data over an image space, converting the cone-beam data to parallel fan-beam data or 
the fan-beam data to parallel-beam data with respect to a rotation angle v^thin the image space, 
performing a shift variant filtration of the parallel-beam data or parallel fan-beam data within the 
image space and converting the processed data to images through backprojection or other means. 



Brief Description of the Drawings 

[0006] FIG. 1 depicts a tomographic image reconstruction system in accordance with 
an illustrated embodiment of the invention; 

[0007] FIG. 2 depicts a flow chart of process steps that may be used by the system of 

FIG. 1; 

[0008] FIG. 3 depicts fan-beam geometry that may be used by the system of FIG. 1 ; 

[0009] FIG. 4 depicts images reconstructed from fan-beam sinograms; 

[0010] FIG. 5 depicts contours of the moduli of the Fourier transforms that may be 
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used by the system of FIG. 1; 

[0011] FIG. 6 depicts average moduli of the Fourier transforms of reconstructed 

images; 

[0012] FIG. 7 depicts variance images comparing the FFBP (first column), the new 
hybrid algorithm (second column) and the previous hybrid algorithm (third column); 

[0013] FIG. 8 depicts theoretical profiles on the center rows in variance images; 

[0014] FIG. 9 depicts variance images obtained for a fan-beam configuration; 

[0015] FIG. 10 depicts images of a bead phantom that may be generated by the 
system of FIG. 1; 

[0016] FIG. 1 1 depicts images of a thorax/shoxilder phantom that may be generated 
by the system of FIG. 1 ; 

[0017] FIG. 12 depicts the images of FIG. 1 1 through a grayscale window 

[0018] FIG. 13 depicts locations at which the resolution properties of two algorithms 
used by the system of FIG. 1 may be evaluated; 

[0019] FIG. 14 depicts AAMs that may be used by the system of FIG. 1 ; 

[0020] FIG. 1 5 depicts contour plots that may be generated by the system of FIG. 1 ; 

[0021] FIG. 1 6 depicts image variances that may be generated by the system of FIG. 

1; 

[0022] FIG. 17 depicts theoretical profiles that may be generated by the system of 

FIG. 1; 

[0023] FIG. 18 depicts a vertically-parallel fan beam data that may be generated by 
the system of FIG. 1; 

[0024] FIG. 19 depicts images that may be generated by the system of FIG. 1 for a 
Shepp-Logan phantom; 

[0025] FIG. 20 depicts comparative profiles that may be generated by the system of 
FIG. 1 using Shepp-Logan information; 

[0026] FIG. 21 depicts AMFTs that may be generated by the system of FIG. 1; 

[0027] FIG. 22 depicts 2D variance images that may be generated by the system of 

FIG.I ; 

[0028] FIG. 23 depicts profiles along the center columns in 2D variance-images that 
may be generated by the system of FIG. 1; 
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[0029] FIG. 24 depicts asymmetric cone-beam configuration information that may be 
used by the system of FIG. 1 ; 

[0030] FIG. 25 depicts images that may be generated by the system of FIG. 1 using a 
Shepp-Logan phantom; 

[0031] FIG. 26 depicts variance-images that may be generated by the system of FIG. 

1; 

[0032] FIG. 27 depicts reconstructed images generated by the system of FIG. 1 using 
asymmetric cone-beam configurations; 

[0033] FIG. 28 depicts fan-beam configuration that may be used by the system of 

FIG. 1; 

[0034] FIG. 29 depicts comparative images that may be generated by the system of 
FIG. l;and 

[0035] FIG. 30 depicts bar pattern profiles that may be generated by the system of 

FIG. 1. 

Detailed Description of an Illustrated Embodiment 

[0036] In this description, mathematical notations are used that characterize the fan- 
beam geometry in a typical CT scanner. As shown in Fig. 3, a continuous fan-beam projection 

q( y, P )\s acquired at a view angle p . The collection ofq(YyP)fox\Y\ ^ Ymzx. P ^ [0.2 ;r ) is 
referred to as the continuous fan-beam sinogram, where y^^ is defined in Fig. 3, The format 
followed herein will be to first sunmiarize the FFBP and the previous hybrid algorithm in their 
continuous and discrete forms. Following the summary of FFBP and previous hybrid algorithm, 
the new hybrid algorithm will be discussed. 

[0037] The well-known FFBP algorithm that reconstructs an image directly fi-om a 
fan-beam sinogram q(y,P) is given by 

f dfi^r dycosr qiy,A .^"^ 



hir.-rl (1) 
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where the constant F denotes the known focal length, 2y^^ is the maximum fan angle, h is the 
spatial domain representation of the ramp filter, and I and are functions of the image-space 
coordinates x and ;/ and the projection angle p , which are given by 



L(x,y,p) = [F'^ ^-x^ ^y'^ +2F{xsmP-y(^osfi and 



y;(x,j;,)ff) = arctan 



XCOSyg + ysinyg 



(2) 



F^xsmfi-ycosfi 



It has been demonstrated that the factor L can significantly amplify noise and aliasing artifacts m 
images reconstructed from data acquired with configurations of small focal lengths and/or a large 
maximum fan angle. 

[0038] One prior art hybrid algorithm may be described in a continuous form as 
follows. Let denote a contmuous parallel-beam sinogram, which, in the absence of data 

noise, is identical to the corresponding fan-beam sinogram provided that 



This observation implies that the parallel-beam sinogram can be obtained fi-om the corresponding 
fan-beam sinogram. Let QSy) denote the Fourier series expansion of the fan-beam sinogram 
q{y,P) with respect to p . The previous hybrid algorithm first calculates the parallel-beam 
sinogram fi*om g^O^) (i.e., fi-om the fan-beam sinogram qiy^P))^ 



^ = Fsm/ and ^ = >ff+;^. 



(3) 



= 1 f}e-^-^QMH'^re''"'Q.('r) V 



(4) 
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and then reconstructs an image from the computed parallel-beam sinogram by use of the parallel- 
beam filtered backprojection algorithm as 



(5) 



where 



^0 j^, ^) = X cos ^ + >^ sin ^. 



(6) 



[0039] The new hybrid algorithm will be considered next. FIG. 1 depicts a system 10 
that may be used with the new hybrid algorithm in accordance with an illustrated embodiment of 
the invention. As shown in Eq. (18) below, the previous hybrid algorithm invokes a ID 
interpolation along ;^for estimating the discrete parallel-beam sinogram on uniform grids 
along 4^ . Such a ID interpolation may result in reduced image resolution in situations where fan- 
beam samples along ;^are sparse. Below, a new hybrid algorithm is described that can avoid 
such a ID explicit interpolation and thus retain the resolution property of the FFBP algorithm. 

[0040] Using Eqs. (3) and (6), one can define the quantity y^ipc^y.^') ^ follows 



/o(^'J^^^) = arcsin 



(7) 



Substituting this definition and q\y^^) into Eq. (5) yields 




(8) 



where q\y,<l>) is given in Eq. (4), and 
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h(Fsin/Q-Fsin/)^ \ dv\v\€ 



Jlm^isinrQ-smr) 



(9) 



We now define a new frequency variable as v' = vF(sin^, - sin -y) . Using v' in Eq. (9) 
yields 



;t(F sin ^0-''^ sin y)= 



n-r 



sin/o - sin/ 



n-r 

sin/j, -sin/ 



-]2 

£^'|v'| 



COS 



COS 



^(/o-y)- (10) 



Substituting Eq. (10) into Eq. (8), we obtain the new hybrid algorithm as 



a(^.>') = :^ r '^^f"" dy cos jq'(y,(p) 



n2 



Kn-r)g(royr% (H) 



where 



cos 



COS 



2 J 



(12) 



Inspection of Eqs. (1) and (11) indicates that the function g(;^o,y) changes the shift-invariant 
filtration in the FFBP algorithm into a shift-variant filtration in the new hybrid algorithm (see the 
integrations over yin Eqs. (1) and (11)). However, because the new hybrid algorithm does not 
involve a factor similar to the L factor (see Eq. (2)) in the FFBP algorithm, it is hypothsized that 
the new hybrid algorithm amplifies noise and aliasing artifacts less significantly than does the 
FFBP algorithm when the focal length F is comparable to or smaller than the FOV size. This 
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hypothesis is quantitatively verified in numerical studies below. 

[0041] Equations 7-12 provide the basis of the new and novel algorithm used by the 
system 10 of FIG. 1. FIG. 2 is a flow chart of process steps that may be used by the system 10 
of FIG. 1. 

[0042] In general, a sampling system 12 may collect fan-beam samples 
conventionally. A menu graphical user interface 16 may be used to enter values and to provide 
the overall control of the system 10, A display 18 may be used to display the images provided 
by the system 10. 

[0043] Once collected, a fast Fourier processor 20 performs a fast Fourier transform 
on the collected data within a Fourier processor 20. Once transformed, a linear combination of 
complementary data elements may be formed from the transformed data within a combination 
processor 22 as described in U.S. Pat. No. 6,115,466 by the inventor of the instant application 
(incorporated herein by reference). 

[0044] Follov^ng formation of the linear combination of data elements, the linear 
combination may be filtered in the spatial domain by a spatial filter 26. Filtering in that spatial 
domain may be accomplished with the use of a shift variant filter 28. The variation in the 
filtering function may be accomplished using a squared trigonometric function 30 as shown in 
Eq. 12 in combination with the arithmetic processor. 

[0045] The reconstruction algorithms in discrete forms may be considered next. The 
FFBP, the previous hybrid, and the new hybrid algorithms in their continuous forms (i.e., Eqs. 
(1), (5), and (11)) are mathematically equivalent, but only in their continuous forms. It can be 
shown that the noise properties in images reconstructed by use of these algorithms in their 
continuous forms are also identical. In reality, one can reconstruct images by use of algorithms 
only in their discrete forms. Because the FFBP and hybrid algorithms in their discrete forms 
propagate data noise and/or aliasing artifacts differently, they reconstruct images with different 
noise (and/or aliasing) properties. The FFBP and hybrid algorithms may be derived in their 
discrete forms and used to study the noise and aliasing properties in reconstructed images. 
Throughout this description 7+1 and N-^l may be used to denote the number of detector bins at 
a projection view and the number of projection views, respectively. For notational convenience 
in expressing these reconstruction algorithms in their discrete forms, both / > 0 and N> 0 3XQ 
assumed to be even integers. However, similar results can readily be derived for / and N being 
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odd integers. 

[0046] The FFBP algorithm in discrete form may be considered first. Let 
A;. = ^and ;', = /Ar, where i = -//2,-//2 + l,...,//2-l,//2, let A;ff = ^ and p„=n^p, 
where » = 0,1,...,A'^, and let {^(/^y^JIdenote samples ofthe fan-beam sinogram. 

[0047] For a given point (x, y) in the image space and a projection angle fi , one can 
calculate 



^ = int 



Ay 



wv i; ft ^-li^MiPjl k 
Ay 



(13) 



where the function "int [b]" yields the largest integer that is smaller than the number b, and 
r'oix,y,p„)is given by Eq. (2). Therefore, the image aj,(x,>') reconstructed by use of the FFBP 
algorithm in its discrete form can be expressed as 



2 ^ [L(x,y,fi„)] 



■{[l-Tjix,y,P„)]qirM + T](x,yM 



(14) 



where q(y,„p„) denotes the filtered discrete fan-beam sinogram, which is given by 



in 



(15) 



and /f denotes the filtration matrix with elements defined as 
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(16) 



[0048] It may be noted in passing that the two terms in the curly brackets of Eq. (1 5) 
(and in Eqs. (20) and (25) for the previous and new hybrid algorithms discussed below) indicate 
that a linear mterpolation has been used in the baclq)rojection step. One can readily use other 
interpolation schemes to replace the linear interpolation. However, such a change of 
interpolation schemes would result in no fundamental change of the theoretical derivation and 
results discussed herein. 

[00491 The previous hybrid algorithm may be considered next in its discrete form. 

Let A# = ^^f^ and ^, =/A^, where j = -//2,-//2+l,...,//2-l,//2,andlet A^ = Jff and 
^„ = /jA^ where n =0,1,..., AT. Using Eq. (4), one can estimate the discrete parallel-beam 
sinogram p(^j,^„) sampled on uniform grids { }. However, because of the non-linear 
relationship between ^ and (see Eq. (3)), uniform grids { ^, } do not correspond to uniform 
grids } on which the fan-beam sinogram is sampled. For a given ^, we define 



a(^,) = arcsin 
/ = int 



'^1 



Ay 
Ay 



(17) 



In the previous hybrid algorithm, the discrete parallel-beam sinogram ^(^, , ^„ ) on uniform grids 
{ } is first estimated from the discrete fan-beam sinogram sampled on uniform grids {p', } by 
invoking a ID linear interpolation as 

p{^M = {^-K^MiYM^K^y^iY,.xA\ (18) 
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where, according to Eq. (4), 



1 Nil 
^ iif=-A^/2 

1 Nil 

q'iyuM=): Z [e-^'"''^''^e:(r/.,)+(-i)"e^'"''<^'^e:(-r,.,)y'"'", (19) 



and denotes the discrete Fourier transform obtained from the discrete fan-beam sinogram 

q{Yi , ) with respect to . Therefore, the image ap„ (x,y) reconstructed by use of the previous 
hybrid algorithm in its discrete form can be expressed as 

where 



^ = int 



^4 



(21) 



^0 (.X' y> ) is given by Eq. (6), and p(^^ , ^„ ) denotes the filtered discrete parallel-beam 
sinogram, which is obtained from the interpolated discrete parallel-beam sinogram in Eq. (19) as 

pi4M = X m><i>„m. -4,)- (22) 

/=-//2 
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[0050] The new hybrid algorithm may be considered next in its discrete form. Let 
Ay = %^and y,. =/Aj^,where i = -//2,"//2 + l,...,//2-l,//2,andlet = ^and<i>„ =«A^, 
where w= 0,1, ...,N. For a given point (x,y) in the unage space and a projection angle one 
can calculate 



^ = int 



Ay 
Ay 



(23) 



where ro(^,y,(^„) is given by Eq. (7). 

[0051] Based upon Eq. (4), a discrete version of the quantity q'(y,(p) can be obtained 

as 



1 N/2 



Considering Eqs. (1 1) and (24), the image a^^^j (x^y) reconstructed by use of the new hybrid 
algorithm in its discrete form can be expressed as 



A^ ^-1 _ _ 

^NH{^^y)^^Y^{{^- K^^yA)]<iXrM^ (25\ 

11=0 ^ ^ 



where 

111 

q'(rM = 2^osr>q'(rMH(k,i)G(k,i), (26) 

iW/2 
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H(k, i) is given by Eq. (16), and G(k i), which is the discrete form of the function gi/o/) in Eq. 
(12), is given by 



cos 



COS 



l2 



(27) 



2 J 



[0052] The noise properties of the reconstructed images will be considered next. 
In the presence of noise, the acquired fan-beam sinogram q(/fy/3„)can be treated as a stochastic 
process. Throughout this description, a bold letter and the corresponding normal letter is used to 
denote a stochastic process and its mean, respectively. It has been shown that, to a very good 
approximation, data noise in CT can be treated as an uncorrelated stationary stochastic process. 
In this situation, the covariances of the measured fan-beam sinogram q(/,.,)5„) can be expressed 
as 

cov{q{r,,fiMiyr>A)} = 9o^(r/ -rrWn - AO. (28) 



where S (.) denotes the Kronecker delta function, and qo is a constant denoting the noise level. 
Consequently, images reconstructed from and quantities estimated from q(y,.,yff^.) can also be 

treated as stochastic processes. Below, we derive image variances obtained with the FFBP, the 
previous hybrid, and the new hybrid algorithms in their discrete forms. 

[0053] The noise properties of the FFBP algorithm in discrete form may be 
considered. Using q(y,.,)3^) in Eq. (14), one can reconstruct a stochastic imagea^ (jc,>'). 

Considering Eqs. (14) and (28), one can show that the variance of a^r (x, y) can be expressed as 

N p2 1/2 

t:i4[L(x,y,p„)] 
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X {[1 - rj{x,y,P„)]Hik,i) + T](x,y,J3„)Hik + 



(29) 



where matrix //is given by Eq. (16), and A: and Tj(x,y,fi„)aTe given by Eq. (13). 

[0054] The noise properties of the previous hybrid algorithm in discrete form may be 
considered. Because q (y, , )5„ ) is a stochastic process, its discrete Fourier transform Q'„ (X/ ) with 
respect to fi„ is also a stochastic process. So is the quantity q'(r, .(^») that is calculated from 
Q^(y.)by use of Eq. (19). It can be shown that q'(r<.<*n) is an uncorrected stationary stochastic 
process with 

The result of Eq. (30) indicates that q'(r„A) is an uncorrelated stationary stochastic process. 
Consequently, the discrete parallel-beam sinogram p(^4,^„) estimated from q'(yi,^„) by use 
of Eq. (18) should also be interpreted as a stochastic process. So is the filtered version 
p(^t,^„)ofthe interpolated discrete parallel-beam sinogram p(^t,^„) (see Eq. (22)). 
Therefore, using Eqs. (18), (22), and (30), one can express the covariance ,^„) as 

III in 

Cov{p(4,^„),P(4.«^„')} = ^o(A^)' Z Z^i^i^-^M^'^-l^) 

+Xi4,)(L-H4,mr, -^n-l (31) 

where A(^,)and 1 are given by Eq. (17). 

[0055] Using p (^^ , ^„ ) in Eq. (20), one can reconstruct a stochastic image 
»PH (x, y) . Considering Eqs. (20) and (3 1 ), one can show that the variance of a p„ (x, y) can be 
expressed as 
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Var{ ^^(x,y)} = ^-^ti\.^-ax,y,^n)fV^mM} 

4 /i=0 

+ 2f(x,;;,^)„)[l-^(x,>',^)„)]Cov{p(^„^>jJ(^,,„(!>„)}}, (32) 

where^and C(x,y,^„) aregivenby Eq. (21), Far{p(^t,^„)}is the variance o/ {p(^»,^„) that 
can be obtained by letting A: = A:' and n = «' in Eq. (3 1 ), and Cov{p(^^ , <p„ ), pC^*^, , ^„ )} can be 

obtained by letting = A: + 1 in Eq. (3 1). 

[0056] The noise properties of the new hybrid algorithm in discrete form may be 
considered next. Because QUr;)is a stochastic process, the quantity q'(;',,^„)that is calculated 
from Q'„(/, )by use of Eq. (24) is also a stochastic process. Again, it can be shown [18] that 
q'(y„^„) is also an uncorrelated stationary stochastic process with 

Cov{q\r<J„),q'(r, ,(/>„•)} = q,8iy, -y,)S{(t>„ -(!>„.). (33) 

Considering Eqs. (25) and (33), one can show that the variance of the reconstructed stochastic 
image 2if^„{x,y) can be expressed as 

Var{2i,„{x,y)} = -^m\^<l>fjL f.^^os' {[\ - ti{x,y,(l>„)]H{k,i)G{k,i) 

4F 11=0 i=-// 2 

+ H{x,y,(l>„)H{k + U)Gik + 1,/)}', (34) 



where matrices and G are given by Eqs. (16) and (27), respectively, and iand /i(x,y,^j)aie 
given by Eq. (23). 

[0057] Quantitative results may be considered next. Computer-simulation studies 
were conducted to evaluate the FFBP, the previous hybrid, and the new hybrid algorithms and to 
validate quantitatively the theoretically predicted image variances in Eqs. (29), (32), and (34). 
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[0058] The data was obtained as follows. In the numerical studies, the Shepp-Logan 
phantom on a 128 x 128 array of square pixels was used to generate noiseless fan-beam 
sinograms that simulate data acquired by use of different scanning configurations. A fan-beam 
sinogram is composed of 180 projection views, each of which contains 257 detector bins. Using 
such noiseless sinograms as the means, noisy sinograms were generated by adding uncorrelated 
stationary Gaussian noise. As shown in Fig. 3, a fan-beam scanning configuration can be 
specified by the focal length F and the center fan angle C. In simulation studies, two fan-beam 
configurations were considered, one of which is specified by (a) F = 363 pixels and C = 0.43 jt 
and the other by (b) F = 107 pixels and C = 0.90 ;r . Configuration (a) represents one that is 
used in a typical CT scanner, whereas configuration (b) represents one that has the same FOV as 
that of configuration (a), but has a focal length shorter than that of configuration (a). We first 
reconstruct images from simulated noiseless and noisy fan-beam sinograms of the Shepp-Logan 
phantom. 

[0059] In an attempt to investigate the resolution properties of the FFBP, the previous 
hybrid, and the new hybrid algorithms, the moduli of Fourier transforms of images reconstructed 
by use of the three algorithms may be calculated and compared. Such moduli can provide 
infoimation about the frequency contents of and possible distortions in reconstructed images. 

[0060) The theoretical predictions for image variances reconstructed by use of the 
three algorithms in their discrete forms are given by Eqs. (29), (32), and (34). In an attempt to 
verify the validity of these analytic formulas, empirical variances were calculated from images 
reconstructed from a large number of statistically independent sets of fan-beam sinograms by use 
of these algorithms. Specifically, for each of the two configurations (i.e., configurations (a) and 
(b)), 4000 sets of fan-beam smograms were generated that contain uncorrelated stationary 
Gaussian noise. Using each of the three algorithms, 4000 noisy images have been reconstructed, 
which can be used to calculate the empirical image variances as 



V{x,y) = 



1 

M-l 



M 1 

t3 M 



-i2 



(35) 



.m=l 
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where a„(x, 3/) denotes the mth image in the total of M(= 4000) noisy images reconstructed by 
use of each of the three algorithms. 

[0061] The reconstructed images may be considered next. From fan-beam sinograms 
of the Shepp-Logan phantom generated with configurations (a) and (b), images may be 
reconstructed and displayed, as shown in Figs. 4a and 4b, respectively. Images in the first and 
second rows were reconstructed from noiseless and noisy sinograms by use of the FFBP (first 
colunrn), the new hybrid (second column), and the previous hybrid (third column) algorithms, 
respectively. As shown in Fig. 4a, for configuration (a) with a large F and a small C, images 
reconstructed by use of the three algorithms appear to be visually similar. However, as shown in 
Fig. 4b, for configuration (b) with a small F and a large C, the FFBP algorithm amplifies data 
noise and aliasing artifacts more significantly than do the hybrid algorithms. Images obtained 
with the new hybrid algorithm appear to have a better resolution than do images obtained with 
the previous hybrid algorithm because the former can avoid the ID linear interpolation (see Eq. 
(18)) that is used in the latter. Results in Fig. 4 also suggest that noise levels in images 
reconstructed by use of the new hybrid algorithm is more uniform than that in images 
reconstructed by use of the FFBP and the previous hybrid algorithms. This observation will be 
verified quantitatively below. 

[0062] Image resolution will be considered next. In an attempt to compare the 
resolution properties in images reconstructed by use of the three algorithms, we consider images 
that contain only point-like structures. Fan-beam sinograms were generated with configurations 
(a) and (b) from three images each of which contains a single point on pixel (68, 68), (80, 48), 
and (90, 90), respectively. 

[0063] In Fig. 5, the contours of the moduli of the 2D Fourier transforms of the 
reconstructed points are shown as fimctions of spatial fi'equencies. The first and second rows in 
Fig. 5 display the results obtained for configurations (a) and (b), respectively. The first, second, 
and third colunms in Fig. 5 display the results obtained for points on pixels (68, 68), (80, 48), 
and (90, 90) by use of the FFBP (dotted curve), the previous hybrid (dash-dotted), and the new 
hybrid (solid curve) algorithms, respectively. As shown in Fig. 5, the solid and dotted curves 
virtually coincide with each other, suggesting that, in the absence of data noise, images recon- 
structed by use of the new hybrid and the FFBP algorithms have virtually identical frequency 
contents. It can also be observed that the dash-dotted curves are generally inside the corre- 
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spending solid curve, suggesting that the image resolution obtained with the previous hybrid 
algorithm is inferior to that obtained with the new hybrid algorithm. One can also convert the 
moduli on 2D Cartesian grids in Fig. 5 onto polar grids and then calculate the average moduli 
over the polar angles. Such average moduli (AMs) play a role similar to that of the conventional 
modulation transfer function. Fig. 6 displays such average moduli, which clearly demonstrate 
that the new hybrid algorithm has a resolution property virtually identical to that of the FFBP 
algorithm and better than that of the previously hybrid algorithm. 

[0064] Image variance will be considered next. For each of the two fan-beam 
configurations, 4000 sets of noisy sinograms were generated from which noisy images were 
reconstructed by use of each of the three algorithms. Using these reconstructed images in Eq. 
(35), empirical image variances were computed. The corresponding image variances were also 
calculated by use of the theoretical formulas in Eqs. (29), (32), and (34) The results displayed in 
Figs. 7a and 7b are for configurations (a) and (b), respectively. The first and second rows in Figs. 
7a and 7b display the theoretical and empirical variances obtained by use of the FFBP (first 
column), the new hybrid (second column), and the previous hybrid (third column) algorithms. In 
Fig. 8, the profiles of these variance-images are shown on horizontal rows through the image- 
array centers, which were obtained with the FFBP (dotted), the new hybrid (solid), and the 
previous hybrid (dash-dotted) algorithms for configumtions (a) and (b). The corresponding 
empirical results obtained with the FFBP (+), the new hybrid (*), and the previous hybrid (A) 
algorithms are also displayed. The empirical results in Figs. 8 and 9 clearly confirm our 
theoretical predictions (see Eqs. (29), (32), and (34) for variances in images reconstructed by use 
of the three algorithms. 

[0065] For configuration (a) with a large F and a small C as shown in Figs. 7a and 
8a, variances in images reconstructed by use of the FFBP algorithm appear to be generally 
imiform, except that variances in peripheral regions are slightly higher than those in the inner 
regions in the image space. On the other hand, variances in images reconstructed by use of the 
new hybrid algorithm are virtually uniform throughout the image space. Variances in images 
reconstructed by use of the previous hybrid algorithm are, in general, lower than those obtained 
with the FFBP and the new hybrid algorithms, reflecting the impact of the ID linear interpolation 
in Eq. (18) that converts a discrete fan-beam sinogram to a discrete parallel-beam sinogram. As 
shown in Figs. 5 and 6, the reduction of image variances obtained with the previous hybrid 
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algorithm can, however, result in a reduction of image resolution when samples of the fan-beam 
sinogram are sparse. 

[0066] For configuration (b) with a small F and a large C, as shown in Figs. 7b and 
8b, image variances obtained with the FFBP algorithm are highly non-uniform, because, in this 
situation, the factor I in Eq. (2) can significantly amplify noise and aliasing artifacts in 
peripheral regions in the image space. In contrast, variances in images reconstructed by use of 
the new hybrid algorithm remain relatively uniform throughout the image space. Again, 
variances in images reconstructed by use of the previous hybrid algorithm appear to be non- 
uniform and are, in general, lower than those obtained with the FFBP and the new hybrid 
algorithms because of the additional ID linear interpolation in Eq. (18). However, as shown in 
Figs. 5 and 6, such an interpolation can lead to a reduced image resolution when samples of the 
fan-beam sinogram are sparse. 

[0067] It is also interesting to observe in Figs. 7 and 8 that the level and uniformity of 
image variances obtained with the new hybrid algorithm remain relatively unchanged for the two 
configurations that have distinctly different focal lengths and center fan angles. However, the 
level of non-uniformity in variance-images obtained with the FFBP and the previous hybrid 
algorithms are strongly dependent upon the configuration parameters (e.g., the focal length F). 
For example, as shown in Fig. 7, variance-images obtained with the previous hybrid algorithm 
contain prominent ring patterns. For configuration (a), because F is relatively large (as compared 
to the FOV size), the relationship between ^ and is close to being linear. Therefore, the linear 
interpolation in Eq. (18) results in an overall decrease of image variances from the center to the 
outer regions in llie image space. For configuration (b), however, because F is comparable to the 
FOV size, the non-linearity between ^ and becomes significant. In this case, the linear 
interpolation in Eq. (18) leads to oscillating ring structures within the variance-images in Fig. 7. 

[0068] When the focal length F approaches infinity, the fan-beam configuration 
approaches the parallel-beam configuration. In this situation, the three algorithms (even in their 
discrete forms) should yield images with identical properties. Sinograms were generated by use 
of a fan-beam configuration with a focal length F(= 50000 pixels) that is much larger than the 
FOV size. Therefore, this configuration can, in a practical sense, be regarded as a parallel-beam 
configuration. From 4000 sets of noisy sinograms generated with this configuration firom the 
Shepp-Logan phantom, images were reconstructed and the empirical image variances were 
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computed from the reconstructed images. For this configuration, Eqs. (29), (32), and (34) were 
used to calculate the corresponding theoretical image variances. Both theoretical and empirical 
results of this study are shown in Fig. 9. Again, the theoretical predictions (the first row in Fig. 9) 
are clearly confirmed by the empirical results (the second row in Fig. 9). Also, as expected, 
identical image variances were obtained with the FFBP (first colunm), the new hybrid (second 
colunm), and the previous hybrid (third colunm) algorithms for this configuration. The fine ring 
structures that appear in Fig. 9 can be attributed to the linear interpolation involved in the 
backprojection steps (see Eqs. (14), (20), and (25)) in all three algorithms. In fact, such fine 
structures also appear in forms modulated by the fan-beam effect in the variance-images in Fig. 
7. For example, for the previous hybrid algorithm, as the results in the third column in Fig. 7 
show, such fine structures are visually overwhebned by the more significant ring structures 
introduced by the non-linear relationship between ^ and y and by the linear interpolation in Eq. 
(18). 

[0069] Following is a short discussion of the results of this first embodiment. In this 
work, a new hybrid algorithm for image reconstruction in fan-beam CT has been described. This 
new algorithm is mathematically equivalent in its continuous form to the FFBP and the previous 
hybrid algorithms. However, under realistic discrete conditions, the new and existing algorithms 
respond differently to data noise and aliasing artifacts. Using discrete forms of these algorithms, 
analytic expressions were derived for noise properties (e.g., image variances) in reconstructed 
images. Computer-simulation studies were performed to evaluate and compare the resolution and 
noise properties of the existing and the new algorithms. Both theoretical and numerical results 
clearly indicate that the new algorithm combines the favorable resolution property of the FFBP 
algorithm and the favorable noise property of the previous hybrid algorithm while eliminating 
their shortcomings. 

[0070] As the theoretical and numerical results demonstrated, one of the highly 
desirable characteristics of the new hybrid algorithm is that its noise properties are much less 
dependent upon the scanning configurations than are those of the FFBP and the previous hybrid 
algorithms. For a given FOV, as the focal length decreases, the noise properties of the new 
hybrid algorithm remain relatively unchanged, whereas the FFBP algorithm amplifies data noise 
and aliasing artifacts more significantly, and the previous hybrid algorithm yields variance- 
images with more conspicuous ring structures. These results may have implications for the 
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design of compact CT systems with reduced focal lengths without sacrificing the FOV sizes. For 
a compact CT system with a short focal length (e.g., configuration (b)), but an FOV that is the 
same as that of a typical CT system (e*g., configuration (a)), one can use the new hybrid 
algorithm to obtain images with improved noise properties over those obtained with the FFBP 
and the previous hybrid algorithms. Therefore, in terms of noise properties, the new hybrid 
algorithm has an evident advantage over the FFBP and the previous hybrid algorithms when 
applied to reconstructing images firom data acquired with such compact CT systems. On the 
other hand, for a fixed focal length, one can enlarge the FOV size by increasing the maximum 
fan angle. Again, images reconstructed firom data acquired with such a system by use of the new 
hybrid algorithm would have lower noise levels than that in images reconstructed by use of the 
FFBP algorithm and more uniform noise levels than that found in images reconstructed by use of 
the previous hybrid algorithm. 

[0071] As displayed in Figs. 7 and 8, noise levels in images reconstructed by use of 
the new hybrid algorithm are much more xmiform than those in images reconstructed by use of 
the FFBP and the previous hybrid algorithm. Such noise uniformity can be desirable for detec- 
tion/classification tasks that are based upon CT images. Results of the simulation studies validate 
quantitatively the derived theoretical results for image variances in Eqs. (29), (32), and (34), 
which not only can provide insights into the algorithms' precision in estimation tasks but also 
may be used for evaluating their performance by use of statistical/mathematical model observers 
in detection/classification tasks. Performance assessment of a reconstruction algorithm in a 
model-observer study generally require knowledge of image variances and covariances. An 
approximation of such knowledge can be computed empirically fi'om a large number of 
reconstructed images. However, this computation can be time consuming, and, more importantly, 
the empirical resuhs so obtained can have significant statistical variations. In contrast, our 
derived theoretical expressions for image variances and covariances can readily be used in 
model-observer studies, thus avoiding empirical calculation of image variances and covariances 
and potential statistical variations within such empirically computed quantities. 

[0072] Because the new hybrid algorithm calculates a set of shift-variant filtration 
coefficients (see Eqs. (11) and (25)), it is slightly less efficient computationally than is the 
previous hybrid algorithm. However, the new hybrid algorithm remains computationally more 
efficient than the FFBP algorithm. Finally, the new hybrid algorithm can readily be generalized 
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to reconstruction of images from data acquired in the half-scan fan-beam configuration, 
asynmietric fan-beam configuration, fan-beam configurations with equi-distance detectors, and 
fan-beam configurations with spatially variant focal lengths. 

[0073] In another embodiment, half-scan strategy could be used as a means for 
reducing the scanning time and radiation dose delivered to the patient in fan-beam computed 
tomography (CT). Also, in helical CT, the data weighting/interpolation functions are often devised 
based upon half-scan configurations. In half-scan CT, after appropriately normalizing the 
redundancy in data, the fan-beam filtered-backprojection (FFBP) algorithm can be used for 
reconstructing unages. As discussed above, the FFBP algorithm can be susceptible to data noise 
and aliasing artifacts, because it invokes a spatially-variant weighting factor that can amplify data 
noise and sample aliasing. Furthermore, because of the lack of data over the full angular range 
2n in half-scan CT, the FFBP algorithm yields images with highly non-uniform resolution and 
noise properties throughout the image space. Such non-uniformity of spatial resolution and noise 
could considerably impede an image's utility in both estimation and detection/classification tasks. 

[0074J In this embodiment, the algorithm for full scan CT (developed above) has been 
generalized to reconstruct images in half-scan CT. The resolution and noise properties are 
investigated m images reconstructed from half-scan data by use of the proposed new half-scan 
algorithm of this embodiment and half-scan FFBP algorithms. In particular, analytic expressions 
are derived for image variances and numerical studies are performed to compare empirically the 
resolution and noise properties in images reconstructed by use of both algorithms. Quantitative 
results in the numerical studies demonstrate that the proposed algorithm of this embodiment yields 
images with more uniform spatial resolution and lower and more uniform noise levels than does 
the half-scan FFBP algorithm. Extensive simulation studies also verify the validity of the derived 
analytic expressions for image variances. The analytic results on image-noise properties can be 
used in the calculation of model observers that assess image-quality in detection/classification 
tasks. 

[0075] The projection data g(Y,P) of an object fimction a(x^% which are acquired by 
use of an equi-angular fan-beam CT configuration, as shown in FIG. 3, can be interpreted as a 
function of projection angle p and detector index y. In fiill-scan CT, knowledge of ^(y^P) is 
available in the fiill data space W = { lyl <ymaxy 0<P <27i} , where lynua denotes the 

maximum-fan angle. One can divide the full data space ^into four subspaces: A, B, C, and Z), 
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where^ = { lyl <ymax,0<^<2imax-2i} ,B = { lyl <ymax,2ymax-2i<^^n-2f) , 

C={ lyl <y/»ax,7t-2y<P<7t+2yOTax},andD={ lyl <ymax,7t+2Ymox^P<27t} . Itcanbe 

shown that information about q(y,fli) in subspace A is identical to that in subspace C and that 
information about ^(y,p) in subspace B is identical to that in subspace D. This observation leads to 
the development of the half-scan CT in which data g(y,p) are acquired only withm the union space 
of ^, B, and C. Obviously, the half-scan data contain redundant information because information in 
subspaces .4 and C is identical. Such a data redundancy can be normalized into the form shown 
(below) in Eq. (1) of this embodiment, 

^y.P)=>Ky.P)9(7,P). (1) 

where the weighting function w(y,p) equals 1 in subspace 5 and 0 in subspace D and satisfies 
w(y,P)+M'(-y,P+2y+n:)=l in subspace A and C. 

[0077] The conventional half-scan algorithm reconstructs an image directly from the 
weighted projection ^w(Y,P) by use of the FFBP algorithm. Below, the full-scan algorithm 
(discussed above) has been generalized into a form for image reconstruction from ^h<Y,P)- 

[0078] In this regard, let p(^,(|)) denote the parallel-beam projections. It is well known 
that, when (p=^+y and ^ = F sin ^ ' the parallel-beam projection p(^,(p) is mathematically identical 
to the fan-beam projection q(y,^). For a given y, one can apply the Fourier-shift theorem to the 
weighted projections to obtain the equation, 

1 ^' 

where 2j"'^W=— Uw(y>fi)^~^^<^P denotestheFourier series expansion of ^„(y,p) with respect 

to p. One can interpret qm(y,(f) as defining a set of parallel-beam projections. In a discrete case, 
however, the samples of 9w(y,<p) on grids imiformly distributed on the y -axis are spread 
non-uniformly on the ^ -axis because of the non-linear relationship between t, and y. Therefore, the 
conventional parallel-beam filtered backprojection (FBP) algorithm has to be modified so that it 
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can be directly applied to samples of q^v,iy,(?) on uniform grids on the y -axis for image 
reconstruction. 

[0079] Using ^ = F sin y and ? iw<y,<p) in the parallel-beam FBP algorithm, one obtains 
the expression, 



271 Jmax 



a(x,y)^jd^ J rfyFcos y^mtCy,*) Mi^sinyo-Fsin y), (3) 



0 -ymax 



where h{^)= Jrfvjv|e^^"^ denotes the spatial function of the ramp filter \u\ , and yo is determined 

-09 

by the expression 



(j:cos ^+ysm ^^ 
P 



(4) 



Using the expression of ^ = Fsin y , and = F sin / in Eq. (3) of this embodiment, we 
obtain the proposed algorithm for image reconstruction firom half-scan data as 



2k ymax 



0 -ymax 



YO-Y 



sin (yq-y) 



A(Yo-Y)g(YoY)> (5) 



where 



^(Yo.Y) = 



cos- 



Yo-y 



cos- 



YO+Y 



(6) 



[0080] The proposed algorithm in Eq. (5) of this embodiment eliminates the 
undesirable spatially-variant weighting factor in the FFBP algorithm. However, because g(yo,y) is 
not a function of yo-y, the filtration in the proposed algorithm of this embodiment (i.e. the 
integration over y in Eq. (5) of this embodiment)) becomes shift-variant Therefore, the 



24 



conventional discrete fast Fourier transform (DFFT) cannot be used to compute the filtration in Eq. 
(5) of this embodiment. One may, however, use non-uniform DFFT or matrix multiplication 
techniques to calculate the shift-variant filtration. 

[0081] Noise properties will be considered next. The half-scan FFBP algorithm and 
the proposed algorithm (see Eq. (5) of this embodiment) in their continuous forms are 
mathematically equivalent. In practical situations, however, one must use the algorithms in discrete 
forms to reconstruct images because the acquired data are discrete in nature. It is in their discrete 
forms that these algorithms propagate data noise and sample aliasing differently. 

[0082] The proposed half-scan algorithm of this embodiment in discrete form will be 
considered next. Assume that there are 7+1 detectors at a projection view and that there are 
^+1 projection, views uniformly distributed over 2n. For notational convenience, both />0 and 
N>0 are assumed to be even integers. (The results can readily be generalized to situations where 
/and A^are odd mtegers.) Let a v = ^^'"'^ and yri^y, where /=-//2,-//2+l,... J/2- 1,7/2, and let 

Ap = and where «=0,1,...,A^. The expression, q(y.^ p„) , can be used to denote the 

discrete data in half-scan CT. Therefore, ^ ^ _ ^ for Prt^2ymax' The weighted discrete data 
is given by the expression, 

q^^yi. P/i) = HJn M qiJr (7) 

[0083] The discrete form of the reconstruction algorithm depends upon the 
interpolation scheme that is used in the backprojection step. In the derivation below, we assume 
that a linear interpolation is invoked in the backprojection step. One can readily use other 
interpolation schemes to replace the linear interpolation. However, such a change of interpolation 
schemes would not lead to a fundamental change of the results presented in this embodiment. 

[0084] Let and <pi,=n<d ^, where «=0, 1 . ,J^. For a given point (x^) in the 

^ N + l 

image space and a projection angle (p„, one can calculate the values, 
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Ay 



H(jc,y.*n) = — k, 



(8) 



[0085] For a given Yi, one can calculate the discrete Fourier transform, Qj^^fi)* of the 
weighted discrete data 9w(y/, Pn) with respect to Using Qj^\yd> a discrete version of 

9w(y,9) in Eq. (2) of this embodiment can be computed in the form. 



[0086] Substituting the computed qm{ Y/, ^n) into Eq. (5) of this embodiment, one 
obtains the integration over y in Eq. (5) in the fomi, 



1/2 



qy^ri Yjt, <t>Aj) = Ay J] cos y,. ?wK Y/, ^n) H{ k, i) , 
/=: -7/2 



(10) 



where //denotes the filtration matrix with elements defined as 



Yife-Y/ 



sin (Yjt-Y/) 



cos- 



cos- 



2 J 



(11) 



Using Eq. (10) of this embodiment in Eq. (5) of this embodiment and considering a linear 
interpolation in the backprojection step, the proposed half-scan algorithm can be obtained in its 
discrete form, 



a(x,y) = -y X {[^ J^wKYj^, ♦«) +^(^,J^,<|)«) ^w/<yjt + i, W } . (12) 



n = 0 



Using Eq. (12) of this embodiment, image variances are derived below by use of the proposed 
algorithm in its discrete form. 
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[0087] The variances in reconstructed images will be considered next. In the presence 
of noise, the measured data should be treated as a stochastic process. A bold letter and the 
corresponding normal letter may be used to denote a stochastic process and its mean, respectively. 
In CT applications, the fan-beam projections qiy^, P«) are converted from the collected photons on 

the detectors. The numbers of these photons satisfy the Poisson distribution, and it can be shown 
that, to a very good approximation, 9(7,-, P«) can be treated as uncorrected stochastic processes, 



i.e., 



Cov{q(Ypp/i), q(y/s ft^) } = -====== d^^^nn^ (13) 

>1p(YpP«)p(Y/sP«') 



where 8,j denotes the Kronecker delta function, and p(y/, p«) is the mean number of x-rays on 

detector (/». Therefore, considering Eq. (7) of this embodiment, the normalized data ^^^i' is 
a stochastic process, with a covariance given by 



o>(YpP«)®(y,sp«') 

Cov { qvi.(Y/, P/i), qM^(Y/'> M } = -======&ii&nn' ( 14) 

>ip(YpPw)p(Y,sP«') 



Consequently, Qj^\yi)is also a stochastic process because it is the Fourier series expansion 

coefficient of ^^^^M ^ So is the quantity that is calculated from Qj^'Kyd^ By 

considering the defining expression, 



®(Yp M gKy.', P„0 
" >^P(Y/,P/|)P(Y.,P;,') 



one can show that can be interpreted approximately as an uncorrelated stochastic 

process with a covariance given by 



Cov{q^r(yi>^n), qwr(yfAn') } « 7[<J(Y/,Y/') "Yp -Y/') ( 
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where a( y^, y,') is an approximation of G(y^, y^/, (Jn, P^O by assuming that it is a function 
independent of p. 

[0088] Finally, using Eqs. (9), (12), and (16) of this embodiment, the variances in 
images reconstructed by use of the proposed algorithm in its discrete form may be obtained in the 
form, 

Var{si(x,y) } « —(Ay) (A(|>) ^ ^ cos^ y, [a(y,.,y,) +a( -y,-, -y,) ] 

X { [ 1 - ii{x,y, ^n) ] H{ k, i) + piix.y, ^rt) ^( * + 1, 0 }^ (17) 

where k and ii(xj/,(pj) are given by Eq. (8) of this embodiment, and matrix H is given by Eq. (1 1) of 
this embodiment We have previously derived image variances obtained with the FFBP algorithm 
m its discrete form. Because the half-scan algorithm reconstructs an image directly from the 
normalized fan-beam projections in Eq. (1) of this embodiment by use of the FFBP algorithm, the 
expression for image variances obtained with the half-scan FFBP algorithm is essentially identical 
to that for image variances obtained with the fiill-scan FFBP algorithm. Although such an 
expression is not given here, it was used to compute image variances for comparison with the 
empirical knage variances obtained in our numerical studies below. 

[0089] The results of this study will be considered next. Computer-simulation studies 
were performed to evaluate the numerical properties of the proposed and half-scan FFBP 
algorithms. In an attempt to reveal the possible asymmetric resolution and noise properties in 
images obtained with the half-scan reconstruction algorithms, a bead phantom (as defined in Table 
1) and a thorax/shoulder phantom were for assessment of the algorithm performance. (The 
thorax/shoulder phantom corresponds to a 2D slice at z=150 mm within the 3D thorax/shoulder 
phantom that is available at the website: http://www.imp.uni-erlangen.de.) 
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Table 1 



Parameters for the bead phantom, (xoj'o) and radius are in unit of half of the FOM size.] 





(xo^o) 


radius 


intensity 


1 


(0,0) 


0.9 


1.0 


2 


(0,0) 


0.7 


1.2 


3 


(0,0) 


0.08 


2.0 


4-7 


(± 0.566, ± 0.566) 


0.08 


2.0 


8 


(-0.8,0) 


0.08 


2.0 


9 


(0,-0.8) 


0.08 


2.0 


10 


(0.8,0) 


0.08 


2.0 


11 


(0, 0.8) 


0.08 


2.0 


12-15 


(± 0.8, ± 0.8) 


0.08 


2.0 



Both phantoms were used to generate noiseless and noisy fan-beam data. Images with sizes of 
300x 300 mm^ are represented as 512x 512 arrays of square pixels. Image properties not only are 
algorithm-dependent, but also rely upon configuration parameters such as focal lengths and 
maximum-fan angles. In our numerical studies, we generated different fan-beam data by use of two 
fan-beam configurations specified by (a) F=850 mm and Ymflx=0.097i and by (b) F=250 mm and 
Yi«ax=0.347C, respectively. Parameter set (a) specifies a configuration that has large focal length and 
small maxunum-fan angle, whereas parameter set (b) defines a configuration that has a field of 
measurement (FOM) as big as that in configuration (a), but has a much smaller focal length and/or 
larger maximum-fan angle. A fan-beam data set is composed of 1024 projections uniformly 
distributed over 27c, each of which contains 1201 detector bins. Noisy data were obtained by 
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computing the logarithm of the ratio between the entrance fluence and detected fluence. The 
entrance fluence per detector was assumed to be 1.5x10^ The collected fluence by each detector is 
assumed to follow the Poisson statistics with a mean that is the product of the entrance fluence and 
the exponential of the negative fan-beam projection. The redundant information in the half-scan 
data was normalized by use of the weighting function proposed by Parker, which is given by the 
ejq)ression, 



(o(y,p) = < 



„ (•«) 



0 (y,P)eZ) 

Reconstructed images may be considered next. In Figs. 10a and 10b, images are shown of the bead 
phantom reconstructed from half-scan data acquired with configurations (a) and (b), respectively. 
Images in the first and second rows were reconstructed from noiseless and noisy data, respectively, 
by use of the FFBP algorithm (the left column) and the proposed algorithm (the right column). 
For data acquired with configuration (a), as shown in FIG. 10a, both algorithms reconstruct 
noiseless rniages with visually identical quality. But, in the presence of noise, the half-scan FFBP 
algorithm yields a slightly noisier image than does the proposed algorithm of this embodiment. In 
contrast, for configuration (b), the lower-lefl region in the image reconstructed from noiseless data 
by use of the half-scan FFBP algorithm appears to be slightly inferior to that in the image 
reconstructed by use of the proposed algorithm of this embodiment. More prominently, in the 
presence of noise, the half-scan FFBP algorithm yields considerably noisier images than does the 
proposed algorithm of this embodiment. One can see that there are many streak artifacts in both 
images, which are caused by the high attenuation paths along those beads. It is interesting to 
observe that the appearance of these streak artifacts are very different in the images obtained by the 
two algorithms. In particular, in the lower and left regions of the image obtained with the half-scan 
FFBP algorithm, the streak artifacts and data noise have been amplified much more significantly 
than those m the upper and right regions, which makes the structures in the lower and left regions 
very noisy. Conversely, such highly non-uniform amplification of data noise and the streak 
artifacts does not appear in images reconstructed by use of the proposed algorithm. 
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[0090] Fig. 1 la and FIG. 1 lb also display images reconstructed by use of the two 
algorithms from noisy half-scan data of the thorax/shoulder phantom acquired with configurations 
(a) and (b), respectively. Observations and conclusions similar to those for the bead phantom 
discussed above can be made. In an attempt to demonstrate the noise properties of reconstructed 
images, FIG. 12 shows the same images as those in FIG. 11 , but with a much narrower display 
window of [-100HU, 200HU]. 

[0091] Resolution properties may be considered next. Quantitative evaluation for the 
resolution properties of the half-scan FFBP algorithm and the proposed algorithm of this 
embodiment were also conducted. Specifically, phantoms were used that consist of point-like 
structures to generate fan-beam data and reconstruct images from that data. From these 
reconstructed images, the moduli of their Fourier transforms were calculated. To facilitate an 
effective comparison of resolution properties throughout the unage space, a figure of merit may be 
introduced and used to refer to as the areas under the average moduli (AAM). The AAM may be 
obtained by (1) converting the Fourier moduli on the two-dimensional Cartesian grids to polar 
grids, (2) calculating the average modui over the polar angles, and (3) computing the areas under 
these one-dimensional average moduli. In FIG. 13, we display spatial locations on the two 
diagonal lines at which the AAMs are calculated. 

[0092] In Fig. 14, the AAMs obtained by use of the half-scan FFBP algorithm and tiie 
proposed algorithm were plotted as fimctions of locations on the two diagonal lines. It can be 
observed in Figs. 14a that, for configuration (a), images obtained by use of the two algorithms have 
comparable resolution properties except that the resolution property of the half-scan FFBP 
algorithm is slightly less uniform. In contrast, for configuration (b), as Fig. 14b demonstrates, the 
resolution properties of the half-scan FFBP algorithm become highly non-uniform, whereas the 
resolution properties of the proposed algorithm remain virtually uniform. It should be pointed out 
that the seemingly high values of AAM in Figs. 14b obtained with the half-scan FFBP algorithm 
do not imply a high resolution. Instead, such high values can be attributed to the strong artifacts in 
images obtained with the half-scan FFBP algorithm. In Fig. 15, the contour plots of the Fourier 
transform moduli of the reconstructed point-like structures are shown at (I) (-102mm, -102mm) 
and (II) (102mm, 102mm). Clearly, the highly non-circular contours suggest that, for the 
configuration (b), the half-scan FFBP algorithm yields images with more significant distortions 
than does the proposed algorithm. 
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[0093] Noise properties will be considered next. The analytic expression for image 
variances obtained with the proposed algorithm in its discrete form is given by Eq. (1 7) of this 
embodiment. An expression was also derived for image variances obtained by use of the half-scan 
FFBP algorithm in its discrete form, which highly resembles the expression m Eq. (30) in the 
previous embodiment and is thus not given here. Using these analytic expressions, we calculated 
the image variances for configurations (a) and (b). In an attempt to verify the validity of these 
theoretical results, we use a circular water phantom with a radius of 150 mm to generate 4000 
statistically independent sets of noisy fan-beam data. Each set of data is obtained by assuming an 
entrance fluence of 2.5x1 0^ Using the two algorithms, we reconstructed 4000 noisy images, which 
are subsequently used for computing empirically image variances. We display in FIGs. 16a and 
16b theoretically predicted (top row) and empirically computed (bottom row) image variances for 
configurations (a) and (b), respectively. The results obtained with the half-scan FFBP algorithm 
and the proposed algorithm are shown in the left and right columns, respectively. Additionally, in 
FIG. 17, profiles of these image variances are displayed on their central-horizontal rows. These 
results clearly indicate that the theoretical and empirical results agree with each other, confirming 
the validity of the theoretical predictions of image variances. They also demonstrate that the 
proposed algorithm generally yields images with lower and more uniform noise levels than does 
the half-scan FFBP algorithm. The cause of the half-scan FFBP algorithm reconstructing images 
with highly non-uniform noise levels can be attributable to the combination of the spatially-variant 
weighting factor in its backprojection and the lack of data at certain views in a half-scan 
configuration. For example, for the half-scan configurations considered, it can be shown that, for 
projection views around the lower-left region in the image space, the spatially-variant factor is 
generally larger than 1 in the lower-left region and smaller than 1 in the upper-right region of the 
image space, thus resulting in significant amplification and suppression of noise in these regions, 
respectively. This amplification-suppression phenomena can be observed in image-variance 
results, as shown in Figs. 16 and 17. It can be mentioned in passing that the decreasing 
image-variances toward the peripheral regions of the FOV in FIGs. 16 and 17 obtained with the 
proposed algorithm reflect the property of Poisson statistics of the projection data. For Poisson 
data, the variances are determined by their means. For the used water phantom, these means are 
higher m the center portion of the detector array at each projection view, consequently leading to 
the higher variances in the inner region of the image space. We have also performed analytic and 



32 



numerical studies by use of data containing uncorrelated stationary Gaussian noise. These results, 
which are not shown here, demonstrated that, although the half-scan FFBP algorithm still produces 
images with highly non-uniform noise properties, the proposed algorithm yields images with 
uniform noise properties. 

[0094] In half-scan fan-beam CT, the existing FFBP-based algorithm can lead to 
images with highly non-uniform resolution and noise properties in situations where the focal 
lengths are small relative to the size of the FOM and/or the maximum-fan angles are large. In this 
embodiment, an algorithm has been described for image reconstruction in half-scan fan-beam CT. 
Theoretical analysis and numerical studies suggest that the proposed algorithm yields images with 
more uniform resolution property and lower and more uniform noise levels than does the half-scan 
FFBP algorithm. Analytic expressions have been derived for image variances obtained by use of 
the half-scan FFBP and the proposed algorithms in their discrete forms and verified the validity of 
these expressions by use of numerical studies. The proposed algorithm would be particularly useful 
for image reconstruction from data acquired by use of configurations with short focal lengths 
and/or large maximum-fan angles, which may be encountered in compact micro-CT and radiation 
therapeutic CT applications. The derived analytic expressions for the image-noise properties can be 
used for image-quality assessment in detection/classification tasks by use of model-observers. The 
proposed algorithm can also be generalized to other scanning configurations such as fan-beam 
configurations with an ofF-set detector, spatially varying focal lengths, or equispace detectors. The 
strategy discussed here may also be applied to reconstructing images from data acquired with 
cone-beam configurations. 

[0095] In another embodiment, it is know that cone-beam data acquired with a 
circular orbit has insufficient information for exact 3D reconstruction. In an attempt to obtain 
sufficient data sets, non-planar and non-circular acquisition orbits in addition to the circular orbit 
have been investigated. Exact algorithms have been developed for image reconstruction from 
cone-beam data acquired with non-circular and non-circular obits. Despite this, the cone-beam 
scanning configuration with a circular orbit remains one of the most popular configurations and 
has been widely employed for data acquisition in, e.g., micro-CT and in CT imaging in radiation 
therapy because of its minimal mechanical complexity. The Feldkamp-Davis-Kress (FDK) 
algorithm, which is a 3D generalization of the conventional 2D fan-beam filtered backprojection 
(FFBP) algorithm, was developed for reconstructing 3D images from cone-beam data acquired 
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vsdth a circular orbit. Other algorithms such as the T-FDK algorithm have also been proposed for 
improving image quality over those obtained with the FDK algorithm. 

[0096] In the FDK algorithm, the spatially-variant weighting factor can not only 
increase the computational load, but also amplify data noise and aliasing artifacts, especially in 
situations where the focal lengths are comparable to the size of field of measurement (FOM). It 
has been shown that the weighting factor can be eliminated by rebinning the cone-beam data into 
parallel beam data. However, the rebinning mterpolations may lead to a loss of spatial resolution 
when cone-beam data samples are sparse. Through the use of a Fourier rebinning scheme and a 
spatially shift-variant filtration scheme, the algorithm described in previous embodiments can 
eliminate the weighting factor in the FFBP algorithm and avoid an explicit interpolation. 
Nimierical results have demonstrated that the algorithm can improve the noise properties without 
compromising the spatial resolution in reconstructed images. Motivated by these observations 
for the above described 2D algorithm, a circular orbit algorithm is presented in this embodiment 
with spatially shift-variant filtration for image reconstruction fi-om cone-beam data acquired with 
a circular orbit. To a certain extent, the algorithm resembles conceptually the T-FDK algorithm 
in that both entail the conversion of the cone-beam data into parallel fan-beam data and 
reconstruction fi-om such converted data. However, because the new algorithm uses the Fourier 
shift theorem to achieve the data conversion and a spatially shift-variant filtration in 
reconstruction, it can avoid the explicit interpolations involved in the T-FDK algorithm. 
Computer simulation studies have been performed to compare the proposed, FDK and T-FDK 
algorithms. Numerical results of these studies demonstrate that the proposed algorithm has 
resolution comparable to and noise properties that are better than the FDK algorithm. As 
compared to the T-FDK algorithm, the proposed algorithm reconstructs images with an 
improved in-plane resolution. 

(0097] The theory of the FDK algorithm will be considered first. Let / (x,;;, z) 
denote the spatial fimction to be reconstructed, and let q{s^h,p) denote the cone-beam data 
acquired with a 2D flat-panel detector and a circular source orbit of a radius of F, where 
p e [0,2;r] is the projection angle and {s,b) are the coordinates of a rotating coordinate system 
on a virtual detector plane containing the center of rotation and mapped from the real flat-panel 
detector with maximum sizes of 2^^ and 26„ along radial and vertical directions, respectively. 
The FDK algorithm can be expressed as shown below. 



34 




(1) 



where 



F(xcosfi + ysmfi) 
F+xsmfi-ycosfi' 



(2) 



F + xsin /3 - ycos fi' 



(3) 



U(x,y,fi) = 



F + xsinfi-ycosfi 
F 



(4) 



[0098] As discussed above, the spatially- variant weighting factor U can significantly 
amplify data noise and aliasing artifacts when the focal length F is comparable to the size of the 



[00991 The proposed algorithm will be considered next. As shown in FIG. 18, (only 
half space z>0 is shown), the cone-beam data acquired from a circular orbit can be rebinned into 
fan-beam data that are vertically parallel to each other. Similar to the situation in 2D fan-beam 
configuration, the cone-beam rotation angle J3 and the rotation angle of the rebinned vertically- 
parallel fan-beam projections satisfy the equation that follows. 



[0100] Therefore, using the Fourier-shift theorem and Eq. (5) of this embodiment, one 
can obtain the vertically-parallel fan-beam data from the measured circular-orbit cone-beam data 
in the form. 



FOM. 



<^ = /? + arctan(^/F) 



(5) 



(6) 



[0101] where (5,6) denotes the Fourier series expansion of the con-beam data 
q{Syb,fi) with respect to fi . Notice that, for these rebinned data, except for those in the mid- 
plane, no two rays corresponding to the same row of detectors are parallel to each other. Instead, 
if mapping the rebinned data of the same row to the virtual detector plane containing the center 
of rotation, one obtains a curved trajectory. Let u denote the height of this trajectory, which can 
be expressed as a function of s and 6, i.e.. 



[0103] Although p{s,u,^) represents parallel-beam projections, it is not uniformly 
sampled over s. Letting 




(7) 



[0102] Therefore, one can express the data function in the {s,u} space as 





(9) 



one can show that 



(10) 
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Finally, using the results above, the new circular orbit algorithm can be obtained which is given 
as follows. 



0 'S„ 



^ F Y 



(11) 



where 



cost's 



(12) 



accounts for the effect of cone-beam open angle y along the rotation axis, 



Giso,s) = 



s„-s 



2 2 
+ S 



(13) 



F(xcos^ + ysin(^) 
^Jf^ -(xcos^ + j^sin^)^ 



(14) 



and 



z^l - (x cos ^ + y sin (t>f 



^JF^'- {x cos ^-^y sin +xsin^>->'C0S^ 



(15) 



[0104] It can be seen that the original FDK algorithm involves a shift-invariant 
filtration over s, whereas the proposed algorithm invokes a shift-variant filtration because G(Sq,s) 
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(in Eq. (1 3) of this embodiment) is not a function of Sq-s. However, the most important feature 

of the new algorithm is that it eliminates the weighting factor U and thus improves the noise 
properties without sacrificing the in-plane spatial resolution. 

[0105] Furthermore, previous results indicated that image intensities obtained by use 
of the FDK algorithm decrease as one moves away from the mid-plane. As the results below 
show, both the T-FDK algorithm and the algorithm proposed in this embodiment can compensate 
for such an intensity drop. Additionally, if we remove the cone-angle weighting factor cos/ in 
Eq, (1 1) of this embodiment, an even better improvement can be obtained. This algorithm, 
which can be referred to as '*the proposed algorithm without cone-angle weighting", suggests 
that the missing data problem can be compensated for by adjusting the preweighting factor in the 
proposed algorithm. 

[0106] The results will be discussed next. Computer simulation studies were 
performed for comparison and evaluation of the FDK algorithm, the T-FDK algorithm, and the 
proposed algorithm. In these studies, a numerical 3D Shepp-Logan phantom represented by a 

2563 matrix was used. The physical length of each side of the volume was assumed to be 300 
mm. The projection data was assumed to be acquired by use of a flat panel detector rotating on a 

circular-orbit with a focal length F = 455 mm, s^^xb^- 262.5 x 262.5 mm^, and the maximum 

cone angle is thus +/- 30^. The detector panel is subdivided into 257 x 257 bins. The cone-beam 
data were generated over 360 projection views uniformly distributed over Itt . 

[0107] In FIG. 19, images for a 2D slice are displayed at y = -37.5 nmi reconstructed 
by use of the FDK algorithm (upper left), the T-FDK algorithm (upper right), the proposed 
algorithm (lower left), and the proposed algorithm without cone-angle weighting (lower right), 
respectively. All images are shown in a gray-scale window [0,60HU]. The corresponding 
profiles along z direction at the position x = 0 mm are shown in FIG. 20. Because the profiles 
obtained from the T-FDK algorithm and the proposed algorithm are close to each other, only a 
single curve is shown for the two algorithms. From these results, it can be seen that the intensity 
drop with increasing cone-angle in the FDK algorithm can be substantially compensated for by 
the use of the T-FDK algorithm and the proposed algorithm. Additionally, as shown in these 
results, the proposed algorithm without cone-angle weighting appears to outperform the other 
algorithms under study. 
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[0108] In an attempt to evaluate quantitatively the algorithms' resolution properties, 
cone-beam data was first generated from a point-source phantom containing two small spheres 
and then images were reconstructed from the cone-beam data. The two small spheres with radii 
of 0.75 mm, located at (0 nmi, 37.5 mm, 0 mm) and (0 mm, 37.5 mm, 37.5 mm), respectively, in 
the field of view. Using the reconstructed images for slices at z == 0 nmi and z = 37. 5 mm, their 
2D Fourier transformations were calculated and subsequently the 2D moduli of these Fourier 
transforms were calculated. Such moduli provide quantitative information about the frequency 
contents (i.e., resolution properties) in reconstructed images. 

[01091 By normalizing these moduli to their maxima and averaging over polar angles 
in the 2D Fourier space, one can obtain single profiles as a fimction of spatial frequency, which 
is referred to as the average moduli in the Fourier transform (AMFT). (The quantity AMFT 
plays a role similar to that of the modulation transfer function (MTF) that is widely used for 
evaluation of the resolution of a linear shift-invariant system.) The AMFTs shown in FIGs. 21a 
and 21b were computed for the slices at (a) z = 0 nrni and (b) z = 37.5 mm, respectively. The 
dotted, dashed, and solid curves show the results obtained by use of the FDK algorithm, the T- 
FDK algorithm, and the proposed algorithm, respectively. These results clearly demonstrate that 
for both slices, AMFT curves obtained with the FDK and proposed algorithms are virtually 
identical and are higher than that obtained with the T-FDK algorithm, suggesting that the T-FDK 
algorithm has a resolution property inferior to that of the FDK algorithm and that the proposed 
algorithm retains the favorable resolution properties of the FDK algorithm. 

[0110] The image noise properties of the three algorithms were also investigated. 
Using a set of noiseless cone-beam projections as the means, 1000 sets of noisy projections were 
generated by addition of stationary Gaussian noise. From the reconstructed noisy images, 
empirical image variances were calculated. In the first and second rows of FIG. 22, the 2D 
variance-images are displayed at z = 0 mm and z = 37.5 mm, respectively, obtained by use of the 
FDK algorithm (left), the T-FDK algorithm (middle), and the proposed algorithm (right). In 
FIGs. 23a and 23b, the profiles are shown along the center columns of the variance-images in the 
first and second rows of FIG. 22, respectively. It can be seen that, the FDK algorithm amplifies 
data noise in the peripheral region more significantly than do the other two algorithms. As 
discussed above, such a noise amplification can be attributable to the spatially-variant weight 
factor U. The noise level obtained with the proposed algorithm appears to be higher than that 
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obtained with the T-FDK algorithm. This is because the T-FDK algorithm involves additional 
interpolation operations that can reduce image noise levels. However, such a noise reduction is 
obtained at the price of compromising the image resolution. 

[0111] In this embodiment an algorithm has been proposed for image reconstruction 
from cone-beam data acquired over a circular orbit. This algorithm eliminates the spatially- 
variant weighting factor in the FDK algorithm by the use of a Fourier rebiiming scheme and 
spatially shift-variant filtration. It resembles the T-FDK algorithm, but can avoid the 
interpolation operations involved in the latter. Simulation studies have been conducted for 
quantitative evaluation and comparison of the proposed algorithm, the FDK algorithm, and the 
T-FDK algorithm in terms of reconstruction accuracy, spatial resolution, and noise properties. 
Numerical results in these studies clearly demonstrate that the proposed algorithm is more 
accurate than is the FDK algorithm. It yields images with noise levels lower than the noise 
levels in images obtained with the FDK algorithm while retaining the favorable resolution 
properties of the FDK algorithm. The proposed algorithm has a resolution property superior to 
that of the T^FDK algorithm. 

[0112] In another embodiment, asymmetric cone-beam CT can be used for increasing the 
size of the field of measurement (FOM) in micro-CT or in megavoltage CT for patient localization 
in radiation therapy. From the data acquired with such configuration, images can be reconstructed 
by use of the conventional FDK algorithm with data weighted by a smooth weighting fiinction. 
With the increase of the FOM size, however, the FDK algorithm may significantly amplify data 
noise and aliasing artifacts because it involves a spatially-variant weighting factor. In this 
embodiment an asymmetric cone-beam algorithm is proposed to reconstruct images from 
asymmetric cone-beam data. This algorithm eliminates the spatially-variant weighting factor 
involved in the FDK algorithm, and thus reconstructs images with improved noise properties. It 
maintains the image spatial resolution by use of a shift-variant filtration scheme to avoid 
interpolation along the radial direction. Additionally, the proposed algorithm can compensate for 
the intensity-drop effect resulting fi-om the missing data problem in the FDK algorithm. This 
algorithm is particularly robust when applied to asymmetric cone-beam CT with large size of FOM 
and/or small focal angles. 

[0113] Let qa{s,b,P) denote the cone-beam data acquired with a displaced 2-D flat- 
panel detector and a circular source trajectory with a radius of F, where p e [0,2;r) is the 



40 



projection angle, and (s,b) are the coordinates of a rotating coordinate system which represents 
a virtual detector plane containing the center of rotation (COR), as shown in FIG. 24. Suppose 
that the detector is displaced such that the cone-beam data exist only when s € [-5^ , 5„ ] . One 
can weight the acquired data by use of a smooth weighting function a? , which equals 0 for 
s € [s^ ] , 1 for 5 € (5^ , 5„ ] , and satisfies 0){s) + 0){-s) = 1 for 5 € [s^ , 5 J . The weighted 
cone-beam data can be expressed as q^{s,b,P) = q^{s,b,P)Q){s) . In this embodiment, the prior 
art weighting function shown below may be used. 
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[0114] The algorithm described in this embodiment converts the weighted cone-beam 
data to vertically parallel fan-beam data by use of the Fourier-shift theorem, which may be shown 
as: 



(2) 



[0115] where ^ denotes the rotation angle of the rebinned vertially parallel fan-beam 
projections and Q^{s,b) the Fourier series expansion of the weighted data q^(s,b,P) with 
respect to p. A data function p(s,u,^) is calculated as follows: 



( u(F' + s'^ ^ 
7:2 '9 



(3) 



which accomplishes tiie mapping of data on the virtual detector along a straight line. 
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[01 16] Using this result in the FBP formula for parallel-beam and applying the pre- 
weighting factor that accounts for the cone angle y , one can obtain a new asymmetric cone- 
beam algorithm as follows: 

f{x,y,z)= \d<l> \ dscosYp{s,u„,<l)) 

[01 17] To evaluate the algorithm proposed in this embodiment, computer simulations 
may be performed to compare images reconstructed by use of the FDK algorithm and the 
proposed algorithm, both using the weighting function in Eq. (1) of this embodiment. A large 
cone-angle was used to demonstrate the effect of the proposed algorithm and the FDK algorithm. 
The low contrast 3D Shepp-Logan phantom (in a cubic volume with side length of 300 mm) was 
used to generate the projection data with an asymmetric configuration, which has a focal length 
F=283.2 mm, virtual detector size = b„ =337.5 mm with a displaced side length j„=84.4 nun 

(corresponding to maximum cone angles of 50° and -170), 257 x 257 detector bins and 360 
views over Ik . Images are displayed in FIG. 25(a) for 2D slices at y=-37.5 mm (fu"st column), 
z = 0 mm (second column), z=l 8.8 mm (third column), and z=37.5 mm (fourth column) in the 
Shepp-Logan phantom. The first, second and third rows show the corresponding slices in the 
original phantom, the images reconstructed by use of the FDK algorithm and the images 
reconstructed by use of the proposed algorithm, respectively. Display window was set to 
[0,60HU]. FIG. 25(b) also shows the profiles along the center column (z direction) in the first 
column of FIG. 25(a). It can be seen that the intensity drop with the increase of cone-angle in 
the FDK algoithm has been significantly compensated for by the proposed algorithm, therefore, 
the reconstructed area is enlarged along the vertical direction. 

[0118] To investigate the noise properties of the proposed algorithm and the FDK 
algorithm quantitatively, empirical variances were calculated ftora. images reconstructed with the 
two algorithms firom 1000 sets of asymmetric cone-beam data that contain uncorrelated noise. 
The data were generated by use of a sphere water phantom with a radius of 1 50 mm. The noise 
data were added by computing the logarithm of the ratio between the entrance fluence and 
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detected fluence. The entrance fluence per detector was assumed to be 1 .5 x 10^. The detected 
fluence at each detector for each detector for each projection was assumed to follow the Poisson 
statistics for which the mean value was calculated from the exponential attenuation of the 
entrance fluence. In the first and second rows of FIG, 26(a) the 2D variance-images are 
displayed at z = 0 nmi and z = 37.5 mm, respectively, obtained by use of the FDK algorithm 
(left) and the proposed algorithm (right). In FIG. 26(b-c), the profiles are shown along the 
diagonal lines of the variance-images in the first and second rows of FIG. 26(a), respectively. 
One can see that the variances of the images reconstructed by use of the FDK algorithm are 
much higher than those obtained by use of the proposed algorithm, especially in the peripheral 
regions of the FOM. In order to demonstrate the influence of this difference, a bead phantom 
was constructed with structures in the peripheral regions. Asymmetric cone-beam data were 
generated with the configuration described above, noise was added by assuming the entrance 
fluence per detector to be 1 .0 x 10^. FIG. 27 shows the reconstructed images by use of the FDK 
algorithm (upper row) and the proposed algorithm (lower row) for 2D slices at z = 0 mm (left 
column) and z = 37.5 nmi (right column). One can observe that the images reconstructed by use 
of the proposed algorithm is not as noisy as those obtained by use of the FDK algorithm, 
especially in the peripheral regions of the FOM. 

[0119] An asynunetric cone-beam algorithm has been described herein for image 
reconstruction from cone-beam data acquired v^th an asymmetric configuration. This algorithm 
improves the noise properties of the conventional FD algorithm, which amplifies data nose and 
aliasing artifacts significantly when large size of the FOM in the asymmetric configuration is 
intended. A shift-variant filtration scheme is applied under the proposed algorithm to avoid 
interpolation so that the spatial resolution is not sacrificed. The proposed algorithm also 
compensates for the intensity-drop problem in the FDK algorithm, thus increasing the 
application reconstruction regions. 

[0120] In another embodiment, Micro-CT is an important tool in biomedical research 
owing to its capability for noninvasively unaging normal and diseased biology and for monitoring 
disease progress and response to therapy or genetic modifications. Spatial resolution is one of the 
most miportant considerations m micro-CT design. Investigators have proposed to improve the 
image resolution by the use of a small source-to-subject distance (SSD) to increase magnification 
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of the subject under study onto the detector arrays. As a result, improved spatial resolution can be 
achieved without decreasing the size of the detection element. 

[0121] FIG. 28(a) shows a typical scanning configuration with an object completely 
covered by the x-ray beam. One can move the object closer to the focal spot by reducing the SSD, 
as shown in FIG. 28(b). This can lead to a better sampling resolution and consequently a better 
image resolution. In this situation, however, the minimal SSD is determined by the maximum 
fan-angle and the object size. 

[0122] The configuration with asymmetric detector has been proposed for increasing 
the size of field of view (FOV). In this embodiment, we show that it can also be exploited for 
further improvement of the sampling resolution and thus the image resolution. As shown in FIG. 
28(c), a half of the object is now outside of the x-ray beam coverage. In this situation, one can 
observe in FIG. 28(c) that the distance between the center of the FOV (i.e., the object) and the 
focal spot is smaller than the SSD in FIG. 28(b), thus allovsdng further improvement on sampling 
and image resolution. 

[0123] Obviously, one needs to use interpolations to convert the data acquired by use 
of the configuration in FIG. 28(c) into the format so that the conventional fan-beam filtered 
back-projection (FFBP) algorithm can be applied. However, such interpolations can significantly 
compromise the potential gain of image resolution by use of the configuration in FIG. 28(c). 
Additionally, when the FFBP algorithm is applied to data acquired with a configuration with a 
large maximum fan-beam angle, as in FIG. 28(c), it can yield images with considerable noise and 
aliasing artifacts. 

[0124] In this work, instead of using the configuration in FIG. 28(c) for increasing the 
size of FOV, we propose to utilize it for improving the sampling and image resolution. Moreover, a 
new reconstruction algorithm has been developed that can reconstruct images directly fi'om the 
data measured by use of the configuration in FIG. 28(c). Because the resolution enhancement 
algorithm of this embodiment invokes no additional interpolation for data conversion, it can 
produce images with the highest spatial resolution allowable by a fan-beam configuration with a 
fixed maximum fan-angle. 

[0125] The methods used will be discussed first. The term g(s,P) may be used to 
denote the data acquired with the configuration shown in FIG. 28(c) in which the center of rotation 
(COR) (i.e., the center of the FOV) is located on the line connecting the focal point and one end of 
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the actual detector, where Pg [0,27c] is the projection-view angle, and ^ € [ 0, 2s m] denotes positions 

on a virtual detector that passes through the COR and is parallel to the actual detector. The focal 
length F is defined as the distance between the focal spot and the COR. In FIG. 28(c), an axis is 
plotted, denoted ass', across the COR and perpendicular to the line connecting the focal spot and 
the COR. Now the configuration in FIG. 28(c) can be interpreted as a fan-beam configuration that 
has a virtual detector that is parallel to the s 'axis. The relationship between 5 and s is given by 



[0126] where jm is a half of the maximum fan-angle of the x-ray beam. Using this 
relationship, one can convert the data q(s,^) that are measured by use of the actual detector onto the 
grids on the virtual detector (i.e., s ). Subsequently, one can use the FFBP algorithm to reconstruct 
images from the converted data. 

[0127] However, in discrete situations discussed above, such a data-conversion 
involves generally interpolations, which can substantially compromise the image resolution that 
may potentially be gained by use of the asymmetric configuration. In this embodunent, a 
reconstruction algorithm is proposed that can reconstruct images directiy from the measured data, 
thus eliminating the otherwise necessary data interpolations and retaining the gain of sampling 
resolution provided by use of the asymmetric configuration in FIG. 28(c), Mathematically, this 
new algorithm is given by 
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G is defined as 
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and the two Jacobi terais are given by 
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respectively. 

[0128] The results will be considered next. Computer-simulation studies were 
conducted to demonstrate and compare the performance of the FFBP and proposed algorithms on 
image-resolution recovery. In these studies, the three configurations shown in FIG. 28 were 
considered, in which the detector size is 92.5 mm and the distances between the detector and the 
focal spot is 80 mm. However, the focal lengths for configuration (a), (b), and (c) are 45 mm, 
30 mm, and 17 mm, respectively. The FOV sizes are all 30 mm in the three configurations. The 
phantom that were used is composed of several sets of bar patterns with different widths, among 
which the width for the finest bar structures is 0.75 pixel. Three sets of data were generated for the 
three configurations, each of which consists of 600 detector bins and 720 views evenly distributed 
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over 2n. Images reconstructed from these data by use of the FFBP and proposed algorithms are 
displayed in FIG. 29. For a more quantitative comparison, FIG. 30 displays the reconstructed 
profiles of the thinnest bar pattems. As these results show, the proposed algorithm yields unages 
comparable to or better than that obtained with the FFBP algorithm. In particular, for configuration 
(c), it can be observed that the proposed algorithm significantly outperforms the FFBP algorithm. 

[0129] A strategy has been proposed in this embodiment for significantly enhancing 
image resolution without altering the hardware (i.e., the detector, the source, and the relative 
configurations between the detector and source.) In particular, a completely asymmetric fan-beam 
configuration has been described for achieving the highest possible sampling and image resolution. 
Furthermore, an algorithm has been developed for image reconstruction directly from the measured 
data. This algorithm yields images with unproved resolution over that in images obtained with the 
FFBP algorithm and is less susceptible to data noise than is the FFBP algorithm. It should be 
emphasized that the proposed strategy and algorithm can also be applied to cone-beam micro-CT 
system. 

[0130] While a preferred embodiment of the present invention has been illustrated 
and described, it will be understood that changes and modifications may be made therein without 
departing from the invention in its broader aspects. 
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